{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# COVID-19 Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total population, N.\n",
    "N = 1000\n",
    "# Initial number of infected and recovered individuals, I0 and R0.\n",
    "I0, R0 = 1, 0\n",
    "# Everyone else, S0, is susceptible to infection initially.\n",
    "S0 = N - I0 - R0\n",
    "# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).\n",
    "beta, gamma = 0.2, 1./10 \n",
    "# A grid of time points (in days)\n",
    "t = np.linspace(0, 160, 160)\n",
    "\n",
    "# The SIR model differential equations.\n",
    "def deriv(y, t, N, beta, gamma):\n",
    "    S, I, R = y\n",
    "    dSdt = -beta * S * I / N\n",
    "    dIdt = beta * S * I / N - gamma * I\n",
    "    dRdt = gamma * I\n",
    "    return dSdt, dIdt, dRdt\n",
    "\n",
    "# Initial conditions vector\n",
    "y0 = S0, I0, R0\n",
    "# Integrate the SIR equations over the time grid, t.\n",
    "ret = odeint(deriv, y0, t, args=(N, beta, gamma))\n",
    "S, I, R = ret.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXhUVbb4/W8llTkkQEiAUIEQQkLIDISAyIwEBAMy21wRkBa1fbWv3drcX7dD+2irrd3XVvRqK4O2yCAOicwzIiDIEFCQEENCBiAhI5C5qs77xzZFAiEESKUyrM/znKemU6dWDsU6p/bZe22dpmkaQggh2gw7WwcghBCiaUniF0KINkYSvxBCtDGS+IUQoo2RxC+EEG2MJH4hhGhjrJb458+fj4+PD2FhYXW+vmLFCiIiIoiIiOCuu+7i2LFj1gpFCCFEDVZL/HPnzmXTpk03fL1nz57s3r2b48eP89xzz/HII49YKxQhhBA16K214WHDhpGenn7D1++66y7L/UGDBpGVlWWtUIQQQtTQLNr4lyxZwvjx420dhhBCtAlWO+NvqJ07d7JkyRK+++47W4cihBBtgk0T//Hjx1mwYAEbN27Ey8vLlqEIIUSbYbOmnoyMDKZMmcJ//vMfgoKCbBWGEEK0OTprVed84IEH2LVrF3l5eXTu3Jm//vWvVFVVAfDoo4+yYMECvvjiC3r06AGAXq/n0KFD1ghFCCFEDVZL/EIIIZonm1/cFaIlMplMFBQUWH7FCmFLDg4OdOzYEXt7+watL2f8QtyGixcv4uzsjLu7OzqdztbhiDZM0zSuXLlCeXk53t7eDXpPs+jHL0RLU1VVJUlfNAs6nQ53d/db+vUpiV+I2yRJXzQXt/pdlMQvRAv2yiuvEBoaSkREBFFRURw4cMAmcSQlJbFhwwbL48TERF577TVA1e1au3btde/ZtWsXEydObLIYxVVycVeIFmr//v2sW7eOI0eO4OTkRF5eHpWVlTaJJSkpiUOHDnHvvfcCEB8fT3x8vE1iETcnZ/xCtFDnz5+nU6dOODk5AdCpUyd8fX3x9/cnLy8PgEOHDjFixAgAdu/eTVRUFFFRUURHR3P58mUA/v73vxMeHk5kZCSLFi0CIDU1lXHjxtG/f3+GDh3KqVOnAHX2/uijjzJ06FCCgoJYt24dlZWVPP/886xevZqoqChWr17N8uXLeeKJJyyxbtu2rdZ7rlVSUsL8+fOJiYkhOjqahIQEq+03IWf8QtyxF1+0zXbHjh3LSy+9RFBQEGPGjGHmzJkMHz78huu/+eabvPvuuwwZMoQrV67g7OzMxo0b+frrrzlw4ACurq4UFBQA8Mgjj/D+++/Tu3dvDhw4wOOPP86OHTsASE9PZ/fu3aSmpjJy5Eh++eUXXnrpJQ4dOsTixYsBWL58ea3Prus9Nb3yyiuMGjWKpUuXUlRUxMCBAxkzZgxubm63ttNEg0jiF6KFcnd35/Dhw+zZs4edO3cyc+ZMS7t6XYYMGcLTTz/N7NmzmTJlCgaDgW3btjFv3jxcXV0B6NixI1euXGHfvn1Mnz7d8t6KigrL/RkzZmBnZ0fv3r0JCAiw/Bqoz83es2XLFhITE3nzzTcBKC8vJyMjg5CQkFvaJ6JhJPELcYesdcbfEPb29owYMYIRI0YQHh7Oxx9/jF6vx2w2AyqBVlu0aBETJkxgw4YNDBo0iG3btqFp2nU9QsxmM+3btycpKanOz7x2/Yb0KLnZezRN44svviA4OPim2xJ3Ttr4hWihkpOTSUlJsTxOSkqiR48e+Pv7c/jwYQC++OILy+upqamEh4fzpz/9iQEDBnDq1CnGjh3L0qVLKS0tBaCgoAAPDw969uzJ559/DqikXHNq1M8//xyz2UxqaipnzpwhODiYdu3aWa4Z1KWu99QUFxfHO++8Q/V40qNHj97h3hH1kcQvRAt15coVHnroIfr27UtERAQnT57kxRdf5IUXXuCpp55i6NChtYbwv/XWW4SFhREZGYmLiwvjx49n3LhxxMfHM2DAAKKioixNLStWrGDJkiVERkYSGhpa62JrcHAww4cPZ/z48bz//vs4OzszcuRITp48abm4e6263lPTc889R1VVFREREYSFhfHcc89Zaa8JkJINQtyWc+fO4evra+swmtzcuXOZOHEi06ZNs3Uo4hq38p2UM34hhGhj5OKuEKLBru2mKVomOeMXQog2RhK/EEK0MZL4hRCijZHEL4QQbYwkfiFaKHd395uus2fPHkJDQ4mKiqKsrOyWtv/1119z8uRJq8QlbEsSvxCt2IoVK/jjH/9IUlISLi4ut/Te2038ovmTxC9EC7dr1y5GjBjBtGnT6NOnD7Nnz0bTND766CPWrFnDSy+9xOzZswF44403iImJISIighdeeMGyjU8++YSIiAgiIyN58MEH2bdvH4mJiTzzzDNERUWRmpp6w1LNaWlpDB48mJiYGBlx20JIP34h7pSt6jLXcPToUU6cOIGvry9Dhgxh7969LFiwgO+++84y0nbLli2kpKRw8OBBNE0jPj6eb7/9Fi8vL1555RX27t1Lp06dKCgooGPHjsTHx9capTt69Og6SzU/9dRTPPbYY8yZM4d3333XOvtCNCpJ/EK0AgMHDsRgMAAQFRVFeno6d999d611tmzZwpYtW4iOjgZUrZ+UlBSOHTvGtGnT6NSpE6BKM1+rvlLNe/futRSDe/DBB/nTn/7U+H+gaFSS+IW4U7asy/yr6lm4QJVqNhqN162jaRr/8z//w8KFC2s9//bbb9+0tPKtlmoWzZu08QvRRsTFxbF06VKuXLkCQHZ2Nrm5uYwePZo1a9aQn58PYJmFq2ap5fpKNQ8ZMoRVq1YB6mKyaP4k8QvRRowdO5bf/OY3DB48mPDwcKZNm8bly5cJDQ3lz3/+M8OHDycyMpKnn34agFmzZvHGG28QHR1NamrqDUs1/+tf/+Ldd98lJiaG4uJiW/6JooGsVpZ5/vz5rFu3Dh8fH3766afrXtc0jaeeeooNGzbg6urK8uXL6devnzVCEaLRtdWyzKL5ahZlmefOncumTZtu+PrGjRtJSUkhJSWFf//73zz22GPWCkUIIUQNVkv8w4YNq7N3QLWEhATmzJmDTqdj0KBBFBUVcf78eWuFI4QQ4lc269WTnZ2Nn5+f5bHBYCA7O5uuXbs2yvYrKmD1arC3B73++tu6nqv5mqOjWpycat86OIB0YBBCtGQ2S/x1XVpozC5hVVVw5kyjbc5Cp7t6UHB2BlfX2ouLy9X7Hh7Qrp26LwcLIURzYbPEbzAYyMzMtDzOyspq1Itlzs7w4INgMoHR2PDb6qWyUi0VFbVvq6rU/YoK+LWn203p9eoA4OEBnp7QsaNavLzU7S2WUBFCiDtis8QfHx/P4sWLmTVrFgcOHMDT07PRmnlAJdtevRptcxZm89WDQlkZlJZeva25lJSoA8OlS1BeDoWFaqmLqyt06gRdu6qlSxfw9lZNT0II0dislvgfeOABdu3aRV5eHgaDgb/+9a9UVVUB8Oijj3LvvfeyYcMGAgMDcXV1ZdmyZdYKpVHZ2alfE87O6gy+ISorrx4EioqgoADy89VSUKAOFBkZaqlmbw+dO4PBAD16qEWq3Yqa7O3tCQ8Px2g00rNnT/7zn//Qvn17W4d1y0aMGMGbb77JgAEDbul9zz//PMOGDWPMmDG89dZbPPLII7i6ugKqNHT1QLUbSUxM5OTJkyxatOi2Y28s77//Pq6ursyZM4fly5czduxYq3YXtlo/ftEwmqYOCrm5cP48XLigbn8dPFmLl5f6FRMUBP7+6leNsI3m0I+/ZnJ76KGHCAoK4s9//rNNY6pmNBrRN/ALeruJvyZ/f38OHTpkqTfUkMTfXN3u/mgW/fhFw+h06pdDYCAMHQrTp8OTT8KiRTB3LowapZK9o6P6hXDwIHz6Kbz+OqxcCcePq+sNom0bPHgw2dnZlscNLb8McPbsWUaPHk1ERASjR48mIyOD4uJi/P39MZvNAJSWluLn50dVVdUNyzPPnTuXp59+mpEjR/KnP/2JkpIS5s+fT0xMDNHR0ZaRvmVlZcyaNYuIiAhmzpxZ5wQxBw8eZMqUKYDq+u3i4kJlZSXl5eUEBARYPm/t2rW8/fbbnDt3jpEjRzJy5EjLNv785z8TGRnJoEGDyMnJue4zli9fzhNPPGHZ1mOPPcbIkSMJCAhg9+7dzJ8/n5CQEObOnWt5j7u7O3/605/o378/Y8aM4eDBg4wYMYKAgAASExOv2y7AxIkT2bVrl+X9dcX14osv8uabb7J27VoOHTrE7NmziYqKYv369dx///2WbW3dutWyX+6EnDM2U87O6qze3189Npng3DlISYHTp9Uvg+RktTg4qF8BERHQu7dqjhJN58VdL1pnuyMatl2TycT27dt5+OGHAW6p/DLAE088wZw5c3jooYdYunQpTz75JF9//TWRkZHs3r2bkSNH8s033xAXF4eDgwOPPPJIneWZAU6fPs22bduwt7fn//2//8eoUaNYunQpRUVFDBw4kDFjxvDBBx/g6urK8ePHOX78eJ0j9vv168fRo0cBNYtYWFgYP/zwA0ajkdjY2FrrPvnkk/zzn/9k586dljP+kpISBg0axCuvvMKzzz7Lhx9+yF/+8pd692NhYSE7duwgMTGR++67j7179/LRRx8RExNDUlISUVFRlJSUMGLECF5//XXuv/9+/vKXv7B161ZOnjzJQw89RHx8fL2fcbO4pk2bxuLFiy1n/Jqm8Yc//IGLFy/i7e3NsmXLmDdvXr2f0RCS+FsIe3vw81PLqFHqekFyMvz4o7o2cOKEWjw8YMAA6NdPrgm0dmVlZZYSzP379+eee+4Bbr388v79+/nyyy8BVVb52WefBWDmzJmsXr2akSNHsmrVKh5//PF6yzMDTJ8+HftfeyVs2bKFxMRE3nzzTQDKy8vJyMjg22+/5cknnwQgIiKCiIiI6/42vV5PYGAgP//8MwcPHuTpp5/m22+/xWQyMXTo0JvuG0dHRyZOnAhA//792bp1603fc99996HT6QgPD6dz586Eh4cDEBoaSnp6OlFRUTg6OjJu3DgAwsPDcXJywsHBgfDwcNLT0xs9Lp1Ox4MPPsinn37KvHnz2L9/P5988slNP+dmJPG3UB4eEBOjluJi+OknOHJENQft2AG7d0NoqGo+8va2dbStW0PPzBubi4sLSUlJFBcXM3HiRN59912efPLJOyq/DFfH08THx/M///M/FBQUcPjwYUaNGkVJSUm95Znd3Nws9zVN44svviA4OPiGn1GfoUOHsnHjRhwcHBgzZgxz587FZDJZDiT1cXBwsHzGjcpUX6u6tLWdnV2tMtd2dnaW99fcbs31aq6j1+stTWSgDnh3Ete8efO47777cHZ2Zvr06Q2+dlIfaRRoBTw9YcgQeOIJNXahTx/V7fT4cXjvPfj8c3XxWLROnp6evP3227z55ptUVVXdcvnlu+66q1ZZ5eoJXNzd3Rk4cCBPPfUUEydOxN7evt7yzNeKi4vjnXfesQzWrG66GTZsmKV8808//cTx48frfP+wYcN46623GDx4MN7e3uTn53Pq1ClCQ0OvW7dmCWlb8/f3JykpCbPZTGZmJgcPHryl91/7t/j6+uLr68vLL79c63rDnZAz/lZEp1MXgnv1Ut1Gv/sOjh692gwUFgb33KMOFKJ1iY6OJjIyklWrVvHggw/y888/M3jwYEAl8E8//bRW+WV7e3uio6NZvnw5b7/9NvPnz+eNN96wtCNXmzlzJtOnT7dcnAR1cHjsscd4+eWXqaqqYtasWURGRl4X03PPPcfvf/97IiIi0DQNf39/1q1bx2OPPca8efOIiIggKiqKgQMH1vk3xcbGkpOTw7BhwwDVLOTj41Pnr4VHHnmE8ePH07VrV3bu3Hknu/KODRkyhJ49exIeHk5YWNgtVx2eO3cujz76KC4uLuzfvx8XFxdmz57NxYsX6du3b6PEKN05W7niYti7Fw4fVheIHRzg7rvhrrvUfXF7mkN3TtF2PPHEE0RHR1su4NflVr6TkvjbiOJi2LJFnfkDtG8P996regOJWyeJXzSV/v374+bmxtatW2tde7iWJH5xQ+npsHEjVHdrjo6GceNU9VHRcJL4RXMjA7jEDfn7w8KFEBenRv4ePaouAKel2ToyIURTkcTfBtnZweDB6gDg66uagT7+GLZtU72BRMPIj2XRXNzqd1Gaeto4k0n1/tm9WyX9Xr1g6lRVMVTc2MWLF3F2dsbd3b1R55EQ4lZpmsaVK1coLy/Hu4GDdiTxC0C1/X/+uSon3b49zJypSkSLuplMJgoKCiwVZ4WwJQcHBzp27GgZNX0zkviFRXExrFkD2dmq/X/qVAgJsXVUQojGJolf1GI0wvr16qKvTgcTJ0L//raOSgjRmCTxi+toGnz7LVQPgBw1StX8kaZsIVoHSfzihg4dUmf/mgaxsaq/vyR/IVo+SfyiXj//DGvXqt4/MTFqtK8kfyFaNkn84qbOnIHPPlPt/4MHw9ixkvyFaMlkAJe4qYAA1b3T3h7271f1/oUQLZckftEgvXur+YDt7GDPHnXxVwjRMkniFw3Wpw9MmaKaeXbsUF0+hRAtjyR+cUvCwlTffoBvvlHt/0KIlkUSv7hl/furqR7NZli9WqZ1FKKlkcQvbsuYMdC3L1RUqB4/v07vKoRoASTxi9ui08H994PBoOb3XblSdfcUQjR/kvjFbXNwgAceUNU8s7Nh0yZbRySEaAhJ/OKOuLnBjBmqmuehQ5CUZOuIhBA3Y9XEv2nTJoKDgwkMDOS111677vWMjAxGjhxJdHQ0ERERbNiwwZrhCCvx9VWlHADWrYMLF2wbjxCifg0q2ZCbm8vevXs5d+4cLi4uhIWFMWDAAOzsbnzcMJlMBAUFsXXrVgwGAzExMaxcuZK+ffta1nnkkUeIjo7mscce4+TJk9x7772kp6c3yh8mml5Cgurb36EDPPIIuLjYOiIhRF3qPePfuXMncXFxTJgwgY0bN3L+/HlOnjzJyy+/THh4OC+88AKXLl2q870HDx4kMDCQgIAAHB0dmTVrFgkJCbXW0el0lvcXFxc3eIZ40Tzde6+atauwEL7+WlX1FEI0P/r6XtywYQMffvgh3bt3v+41o9HIunXr2Lp1K1OnTr3u9ezsbPz8/CyPDQYDBw4cqLXOiy++yNixY3nnnXcoKSlh27Ztt/t3iGbAwUG193/wASQnw+HDMGCAraMSwno0TcNoNlJlrlK3pirLY5PZhEkzXXdr1sx1vmbWzHU+N9hvMD5uPo0ad72J/4033rjxG/V6Jk+efMPX62pBunZS6pUrVzJ37lz+8Ic/sH//fh588EF++umnepuQRPPWoYMa2bt2LWzeDD16QAPnfxbCqjRNo8pcRWlVKeXGcipNlVQYK6gwVVjuV5oqqTBV1LpfaaqkylRVZ3I3mq3fhznUJ7RpE3+1f/3rX8ybN4927dqxYMECjh49ymuvvcbYsWNv+B6DwUBmZqblcVZW1nVNOUuWLGHTr30ABw8eTHl5OXl5efj4NO4fKZpWWBikpMCxY/DFF7Bgger1I0RjMmtmSqtKuVJ5hSuVV7hccZkrlVcoM5ZRVlVGmbGM0qpSy/2yqjJMmqnR49Db6XGwc1C39g6Wx/Z29tjr7Ou8tdPZ1fmanc7uuue8XRv/zKlB/x2XLl3KU089xebNm7l48SLLli1j3rx59Sb+mJgYUlJSSEtLo1u3bqxatYrPPvus1jrdu3dn+/btzJ07l59//pny8nK85fSwVbj3XsjIUD18duxQNfyFaCiT2cSliksUlRdZlksVl64m+crLlFSWoHFrF5Ic7BxwcXDBWe+Mk70TTnonHO0d673vaO+Ig71Dncldb6e/riWjJWhQ4q9uttmwYQPz5s0jMjKyzqacWhvW61m8eDFxcXGYTCbmz59PaGgozz//PAMGDCA+Pp5//OMf/Pa3v+V///d/0el0LF++vEXuRHE9JyeYOhWWLoV9+yAwUNX1F6Ka0WwkvzSfvNI88krzyC/LtyT5yxWXG5TUXR1ccXd0p51jO9wd3XF3dMfVwRUXBxdc9C64OLiox7/e19vJT09oYHfOefPmkZ2dTVpaGseOHcNkMjFixAgOHz7cFDGKFmz3bjVpu6cnPP64OiCItsVoNpJbksuFKxe4WHLRkuiLyotumNx16PBw8qC9c3s8nT3VrZOnJbm3c2qHm4Mb9nb2TfzXtA4NSvxms5mkpCQCAgJo3749+fn5ZGdnExER0RQxihbMbIaPPoJz51RVz/vus3VEwpqqTFWcv3Kec5fPcf7yeZXsSy9i1szXrWuns6ODcwc6uXaik2snvFy96ODcgfbO7fFw8pCkbkU3TfzFxcVs2rSJ7OxsdDodvr6+xMXF0b59+6aKUbRwubmqi6fJBHPmSJNPa3K54jKZlzLJLM4k81Im5y+fv+4Cqg4dXq5edHXvio+bjyXRd3TpKMndRupN/J988gl//etfGTt2LN26dQNU75ytW7fywgsvMGfOnCYLVLRse/bA9u3S5NPSlVWVkV6UTmphKmcKz1BQVlDrdR06fNx8MHgY6OLeha7tVLJ3tHe0UcSiLvUm/uDgYA4cOHDd2X1hYSGxsbGcPn3a6gGK1qFmk8+AAVdn8RLNm6ZpZF/OJjkvmTOFZzh3+VytdnkneycMHgb8PP3w8/DD4GHASS9H9eau3kvcmqbV2cvGzs7upr16hKjJzg4mT1ZNPocOqUlcpMmneTKZTaQXpXMq7xSn8k5xufKy5TV7nT1+nn4EdAggoEMAvu18sdPJgMuWpt7E/+c//5l+/foxduxYS/mFjIwMtm7dynPPPdckAYrWw8cHhg9X/frXrVNNPjKwq3kwa2bSCtM4nnOc5Pxkyo3lltc8nTzp06kPgR0D6dG+hzTbtAI3vbhbWFjI5s2byc7ORtM0DAYDcXFxdOjQoaliFK2IyQTvvw8XL6qDwMiRto6obcsrzSPpQhLHc45zqeJqwUVvV29CvEPo06kPXd27yviaVqZB3TlzcnJq9erp3LlzU8QmWqmMDDWwy94eHn1Uavk0tUpTJcdzjnP0/FGyL2dbnu/g3IHILpGE+4Tj5eplwwiFtdWb+JOSknj00UcpLi7GYDCgaRpZWVm0b9+e9957j379+jVlrKIV+eYbVb2zRw+YO1fN4Susq6CsgIPZB0m6kGRpynGydyLUJ5TIzpF09+wuZ/ZtRL2JPyoqig8++IDY2Nhaz3///fcsXLiQY8eOWT1A0TqVlcHixVBSApMmQXS0rSNqvTKLM/ku4zuS85Mtz/l5+BHTLYaQTiE42DvYMDphC/Um/t69e5OSklLna4GBgfzyyy9WC0y0fj/+qKp3urjAE0+o+XtF49A0jdTCVPac3cPZ4rOAqiIZ5hNGbLdYurbrauMIhS3V26di/PjxTJgwgTlz5lh69WRmZvLJJ58wbty4JglQtF5hYWpy9tRU2LIF7r/f1hG1fJqmkVaUxvYz2y3t9856ZwZ2G0hst1jcHOXoKhpwcXfjxo0kJCTU6tUTHx/PvdWzawtxBwoK4L33wGiUcg53KrM4kx1pO0grSgPA3dGdwYbBDPAdIIOqRC0N6tUjhDVVl3Po2FH69t+OovIitqRu4eTFkwC46F0Y0n0IA7sNlD73ok71/hcrLi7m1VdfJSEhgdzcXAB8fHyYNGkSixYtkkJtolHcdZdq78/Nhb17Vf9+cXOVpkq+y/iOfZn7MJqNONg5MNhvMHf53YWz3tnW4YlmrN6x1jNmzKBDhw7s3LmT/Px88vPz2blzJ+3bt2f69OlNFaNo5eztYcIEdX/PHtX8I25M0zROXjzJOwfe4duz32I0G4noHMH/F/v/MarnKEn64qZuWqQtOTn5ll8T4nZ89ZWapzcwEGbPlr79dblUcYkNKRs4lXcKAN92vowPHI+fp5+NIxMtSb1NPT169ODvf/87Dz30kGW0bk5ODsuXL7f08hGisdxzDyQnwy+/wKlTEBJi64iaD03TOHTuENvObKPCVIGTvRNjAsYwwHeADLoSt6zepp7Vq1eTn5/P8OHD6dixIx07dmTEiBEUFBSwZs2apopRtBHu7jB6tLq/aRNUVto2nubicsVlPj3+KetT1lNhqqBPpz78buDviOkWI0lf3Bbp1SOalZp1+4cMUb8C2rKTF0/yTfI3lBnLcNG7MDFoIn29+0rCF3fktgtpL1u2rDHjEAJQdfsnTFDt+/v3q54+bVGVqYqEUwmsObGGMmMZgR0DeTzmcUJ9QiXpizt222f83bt3JyMjo7HjEQKA9evhhx/aZhG3/NJ81pxYQ05JDno7PWN7jSXGV5p1ROOp9+JuREREnc9rmkZOTo5VAhICYNQoOHkSzp5Vffxv8FVsdU5ePEnCqQQqTBV4uXgxI3QGnd2lDLpoXPUm/pycHDZv3nzdpCuapnHXXXdZNTDRtrm4qPb9r7+GzZshKAicW3H3dLNmZvuZ7ezN3AtAqHco8cHxUmpBWEW9iX/ixIlcuXKFqKio614bMWKEtWISAoDISDh6VJ3179gBrbU8VIWxgi9+/oLT+aex09kxttdYYrvFStOOsBrp1SOatdxcNVWjpsFvfwu+vraOqHEVlRex8seV5JTk4KJ3YWbYTPzb+9s6LNHK1dur58qVKzfdQEPWEeJ2+fjAoEEq8a9fr7p7thZZl7L48PCH5JTk0Mm1E7/t/1tJ+qJJ1Jv4J02axB/+8Ae+/fZbSkpKLM+fOXOGJUuWEBcXx6ZNm274/k2bNhEcHExgYCCvvfZaneusWbOGvn37Ehoaym9+85vb/DNEazZ8OHh4QHY2HDli62gaR0p+Ch8nfUxJVQm9OvRiQb8FdHTpaOuwRBtx06aeDRs2sGLFCvbu3UthYSF6vZ7g4GAmTJjAww8/TJcuXep8n8lkIigoiK1bt2IwGIiJiWHlypX07dvXsk5KSgozZsxgx44ddOjQgdzcXHx8fBr3LxStwsmTsGZN65it63jOcb4+9TVmzUx0l2gmBk3E3s7e1mGJNuSmlc/vvffe25p05alJ0IkAACAASURBVODBgwQGBhLw68was2bNIiEhoVbi//DDD/nd735n6TUkSV/cSEiIKt72yy+wdStMnmzriG7P91nfs+kX9Sv57u53M7rnaLmIK5rcbY/cvZns7OxahdwMBgPZ2dm11jl9+jSnT59myJAhDBo0qN5mI9G26XQwfrwq4ZyUBC1x7OCu9F2WpB/XK44xAWMk6QubsFrir6sF6dovudFoJCUlhV27drFy5UoWLFhAUVGRtUISLZyXF9x9t7q/bh2YTLaNp6E0TWNn2k52pe9Ch477+9zPYL/Btg5LtGFWS/wGg4HMzEzL46ysLHyv6YtnMBiYNGkSDg4O9OzZk+DgYFJSUqwVkmgF7r4bOnRQ3TwPHrR1NDenaRo70naw++xu7HR2TO07lcgukbYOS7RxN038ZrOZsLCwW95wTEwMKSkppKWlUVlZyapVq4iPj6+1zuTJk9m5cycAeXl5nD592nJNQIi6ODhcHci1cydcumTbeOqjaRrbzmxjT8YelfRDphLmc+v/l4RobDdN/HZ2dkRGRt5yQTa9Xs/ixYuJi4sjJCSEGTNmEBoayvPPP09iYiIAcXFxeHl50bdvX0aOHMkbb7yBl5fX7f0los3o3Vtd7K2sVOUcmqtd6bvYm7kXO50d0/pOI9Qn1NYhCQE0cOTuqFGj+OGHHxg4cCBuNfrRVSdwIZpacTEsXgxVVfDgg9Crl60jqm1f5j62pG7BTmfH9L7TCfGW6cRE89GgxL979+46nx8+fHijByREQ+3dq7p2ennBY4+B/qadk5vG4XOH+eb0NwDc3+d+adMXzU6Da/WcPXuWlJQUxowZQ2lpKSaTiXbt2lk7PiFuyGRSdXwuXlRlnIcNs3VE8FPuT3xx8gs0NMYHjifWEGvrkIS4ToN69Xz44YdMmzaNhQsXAqqP/uSWOoJGtBr29mq2LoBvv4XCQtvGc6bwDF/+/CUaGqN6jpKkL5qtBiX+d999l7179+Lh4QFA7969yW2rc+KJZsXfX03SYjTChg2qmJst5FzJYfVPqzFrZgYbBjO0+1DbBCJEAzQo8Ts5OeHo6Gh5bDQaZcShaDbGjlWTtKSkQHJy03/+pYpLrPhxBRWmCvp692Vsr7Hy/0M0aw1K/MOHD+dvf/sbZWVlbN26lenTp3PfffdZOzYhGsTdXbXxA2zcqLp5NpUKYwUrjq/gUsUlunt2Z0rIFEn6otlr0MVds9nMkiVL2LJlC5qmERcXx4IFC+QLLpoNsxk+/BDOn1f1+8eNa4LP1MysOL6C1MJUvFy8eLjfw7g6uFr/g4W4Qw3u1VNZWcmpU6fQ6XQEBwfXavoRojk4f14lf02DefOge3frft7GlI0cyD6Am4MbC/otoINLh5u/SYhmoEFNPevXr6dXr148+eSTPPHEEwQGBrJx40ZrxybELenaFYYMUYk/IUEN7rKWI+ePcCD7APY6e2aGzZSkL1qUBp3x9+nTh3Xr1hEYGAhAamoqEyZM4NSpU1YPUIhbYTTCBx+ovv1DhsA99zT+Z5wtOssnxz7BpJmYFDyJ6K7Rjf8hQlhRg874fXx8LEkfICAgQCZNEc2SXg+TJqn6/fv2QVZW426/qLyI1SdWY9JMDDIMkqQvWqR6B7l/+eWXAISGhnLvvfcyY8YMdDodn3/+OTExMU0SoBC3ymCAwYNV4k9IgIULG6ecQ5WpilU/raK0qpTAjoGM7TX2zjcqhA3U+9/hm2++sdzv3LmzpWaPt7c3hbYeJilEPUaOVH36L16E7dshLu7Ot7khZQMXrlygo0tHpvWdhp3OatNZCGFVDe7VI0RLk5UFS5eqrp5z5sCdTPVw5PwREpMT0dvpWdBvAV3cuzReoEI0sQYl/rS0NN555x3S09MxGo2W56Uss2judu1Si4eHquDp4nLr2zh/+TxLji7BaDYyuc9korpENXaYQjSpBrV8Tp48mYcffpj77rsPOzv5eStajmHD4Jdf1Nn/+vUwbdqtvb/cWM6aE2swmo3079pfkr5oFRqU+J2dnXnyySetHYsQjc7ODqZMUeWbf/oJgoMhPLxh79U0ja9+/orC8kK6undlfO/x1g1WiCbSoKaezz77jJSUFMaOHYuTk5Pl+X79+lk1OCEay5EjkJgITk6ql0/Hjjd/z3cZ37HtzDac9c4s7L9QBmmJVqNBZ/w//vgj//nPf9ixY4elqUen07Fjxw6rBidEY4mOVtU7f/4ZPv8cHn64/i6e6UXpbD+zHYApIVMk6YtWpUGJ/6uvvuLMmTNSn0e0WDqdGth14YKq6bN589VJXK5VWlVqmUVraPehBHkFNW2wQlhZg67URkZGUlRUZO1YhLAqZ2eYPl3N3PXDD6rN/1qappFwKoHLlZfx8/BjZM+RTR+oEFbWoDP+nJwc+vTpQ0xMTK02funOKVoaX181mGvDBtXm37Wrmqy92qFzh0jOT8ZZ78zUvlNlkJZolRqU+P/6179aOw4hmkxMDJw9CydOwOrVsGABODqq6RM3p24G4L6g+2jv3N7GkQphHTJyV7RJFRWqdn9eHoSEwP1Tq/jo6IfkluTSr2s/4oPjbR2iEFbToMTfrl07y2xblZWVVFVV4ebmxqVLl6weoBDWkp+vkn95OTiFraei0w90cu3EI/0fwdFeOjKI1qtBTT2XL1+u9fjrr7/m4MGDVglIiKbi5aVG8r614hQ//fQD4aH2LBw3VZK+aPVu68rV5MmTpQ+/aBV8/C5h7p0AQFXyPehKuto4IiGsr0GJ/8svv7Qsa9euZdGiRQ2aaH3Tpk0EBwcTGBjIa6+9dsP11q5di06n49ChQw2PXIg7VN1107trGdE9etPZGMuKFSA9l0Vr16Cmnpp1+fV6Pf7+/iQkJNT7HpPJxO9+9zu2bt2KwWAgJiaG+Ph4+vbtW2u9y5cv8/bbbxMbG3sb4Qtx+w6dO0RqYSpujq7898xJfL1GR3o6fPopzJ8Prq62jlAI62hQ4l+2bNktb/jgwYMEBgYS8GsR9FmzZpGQkHBd4n/uued49tlnefPNN2/5M4S4Xfml+WxJ3QLAxKCJtHd1Z9YsWLYMcnJg5UpVw9/BwcaBCmEF9Sb+l1566Yav6XQ6nnvuuRu+np2djZ+fn+WxwWDgwIEDtdY5evQomZmZTJw4URK/aDJmzcxXp76iylxFROcI+nqrkxFnZ5g9G5YsgcxMWLsWZs5UFT6FaE3q/Uq7ubldtwAsWbKE119/vd4N19VLtOZ1AbPZzH//93/zj3/843biFuK27c3YS9alLDycPBgfWLvUsocH/Nd/qQlbkpPh66/VDF5CtCb1nvH/4Q9/sNy/fPky//rXv1i2bBmzZs2q9VpdDAYDmZmZlsdZWVn4+vrW2t5PP/3EiBEjALhw4QLx8fEkJiYyYMCA2/lbhLip85fPszN9JwCTgifh4nD9lFze3vCb38B//gPHj18t8CZn/qK1uOlXuaCggL/85S9ERERgNBo5cuQIr7/+Oj4+PvW+LyYmhpSUFNLS0qisrGTVqlXEx18dDenp6UleXh7p6emkp6czaNAgSfrCqoxmI1+d+gqzZmZgt4H06tjrhuv6+almH0dHOHZM1fWRMe6itag38T/zzDPExMTQrl07fvzxR1588UU6dGhYXXK9Xs/ixYuJi4sjJCSEGTNmEBoayvPPPy/F3YRN7EjbQW5JLl4uXtwTcM9N1+/RQ535OzhAUpIkf9F61Fuywc7ODicnJ/R6fa32eU3T0Ol0UrJBtBhni86yPGk5Op2O+dHzMXgYGvze9HRYsQKqqiAsDCZPrn8SFyGaOynSJlq9CmMF/3fo/ygqL2J4j+G3VWM/PV118ayogIAA1dunRoVyIVoUuVwlWr3NqZspKi+iq3tXhvUYdlvb8PeHuXPB3R3OnIGPP4aSkkYNU4gmI4lftGrJeckcOX8EvZ2e+0Pux97O/ra31bWrGtHboQOcO6f6++flNWKwQjQRSfyi1SqpLCExWXUkGN1zND5u9fdEa4iOHdVE7V27QkGBKut8+vQdb1aIJiWJX7RKmqax7vQ6SqpK8G/vzyDDoEbbtrs7zJsHffuqNv+VK+G776THj2g5JPGLVul4znF+zvsZJ3snJveZ3KBqsrfC0VFN3D5qlEr427apEg/l5Y36MUJYhSR+0eoUlxezIWUDAOMCx1lt7lydDoYNg1mz1IHgxAn44APIyrLKxwnRaCTxi1ZF0zS+PvU1FaYKgr2CieoSZfXP7NMHFi5U7f6FhbB0qTT9iOZN+vGLVuX7rO/Z9Msm3BzceCzmMdwd3Zvss41G2L4d9u9Xj/39IT5eXRAWojmRxC9ajdySXP59+N8YzUZmhc2iT6c+NokjJUVV9SwpUSN8R46EwYOlyJtoPiTxi1bBZDbx0ZGPOH/lPNFdopnUZ5JN4ykthU2bVHVPAF9fuO8+1RwkhK1J4hetwvYz29mTsYcOzh14dMCjOOmbRz2FlBRYtw6Ki9XF4KgoGD1adQkVwlYk8YsWL6M4g2VH1fSg86Ln0d2zu40jqq2iAnbvhgMHwGRSNX6GDoXYWJnaUdiGJH7RolUYK3j/0PsUlhdyd/e7GRMwxtYh3VB+PmzZomb2AmjXTh0A+vWTap+iaUniFy1aYnIiR84foYt7F37b77d3VIunqaSmqgFf58+rxx4eajxAVJQcAETTkMQvWqzkvGRW/rQSe509CwcsbJRaPE1F0+DUKdi1C3Jy1HNubjBwIMTEgKurTcMTrZwkftEilVSW8N4P71FSVUJcrzgG+w22dUi3RdPg5Ek14Kv6F4CDA0REwIAB0gtIWIckftHiaJrG6hOrOZV3ip7tezInck6j1+JpapqmJnvZt0/1BKrWrRv0769m/nJ0tFl4opWRxC9anMPnDvPN6W9wsnfi8ZjH8XT2tHVIjeriRTh0SE3yXl30zcEBgoMhPBwCA8G++V/KEM2YJH7RouSW5PLh4Q+pMlcxJWQKEZ0jbB2S1VRVqWagI0fg7Nmrz7u4qPpAQUHQq5f8EhC3ThK/aDGqTFV8eORDcktyieoSxeQ+k20dUpMpKoKffoIff7x6MRhUL6CePdVBIDhY9RAS4mYk8YsWY93pdRw6dwgvFy8WDliIo33bPNXNzVVjAZKTITu7dhXQTp1UcbjqRUYIi7pI4hctwoncE3x+8nP0dnoW9FtAF/cutg6pWbhyRV0MTk5Wk8BXVtZ+3dsbevRQF4m7dVMHBikWJyTxi2avsKyQDw5/QLmxnHt738vAbgNtHVKzZDKpSeDT09WSkaGuE9Tk6Ki6iHbrpm47dwYvL7lY3NZI4hfNmslsYlnSMrIuZdGnUx9mhs5s8V03m0r1gSAjQ91mZ6trBdeys1O/BHx81NKpk5pDoEMHVVdItD6S+EWztu3MNr7L+A5PJ08eHfAoLg4utg6pRSspUQeAc+fgwgV1vaCw8Mazhbm5qYNAzcXDQy3t2kmJiZZKEr9otk7nn+azHz/DTmfH3Ki5za7qZmtRWQl5eaq3UG4uFBSopbBQzSpWHze3qweC6oOBm5sqOeHmdvW+s7MqSy2aB6sm/k2bNvHUU09hMplYsGABixYtqvX6P//5Tz766CP0ej3e3t4sXbqUHj16WCsc0YIUlBXw78P/ptxYzuieoxnaY6itQ2pzNA0uXVIHgJoHg0uX1HL5MpjNDduWnV3tg4Gzs1qcnBp2Xy5INy6rJX6TyURQUBBbt27FYDAQExPDypUr6du3r2WdnTt3Ehsbi6urK//3f//Hrl27WL16tTXCES1IpamSJUeWkFOSI+36zZjZrJqOqg8ExcXqcfVSWnr1fkXFnX2Wvb0avVxzcXS8/rnq5/V6tdjbq6Xm/fqeu/Z5O7vaS2v5Glqthe7gwYMEBgYSEBAAwKxZs0hISKiV+EeOHGm5P2jQID799FNrhSNaCE3T+Cb5G3JKcvBy8eL+PvdL0m+m7OxU0067dqqXUH2MxusPBOXlarnZ/YoKdaHaZLpawsKWrj0Q1HVwuNnj+pZr1xk6tPGL9Vkt8WdnZ+Pn52d5bDAYOHDgwA3XX7JkCePHj7dWOKKFOJh9kB9zf8TR3pFZYbOazRSK4s7o9VevA9wqTVNJv6rq+qWy8sbPVR8sqhej8dafM5trL1D7flPo16/xt2m1xF9XC9KNztw+/fRTDh06xO7du60VjmgBzhSeYXPqZgAmBU/C283bxhGJ5kCnu9p042LjTl2adjXx17x/o+dutM6tLF2sMFbRaonfYDCQmZlpeZyVlYWvr+91623bto1XXnmF3bt34ySdhtusvNI81pxYg1kzc3f3uwn1CbV1SEJcR6e72v7fklntWnlMTAwpKSmkpaVRWVnJqlWriI+Pr7XO0aNHWbhwIYmJifj4tJzZk0TjKq0q5bMfP6PcWE5IpxBG9xxt65CEaNWslvj1ej2LFy8mLi6OkJAQZsyYQWhoKM8//zyJiYkAPPPMM1y5coXp06cTFRV13YFBtH4ms4k1J9ZQUFZAV/eu3B8iF3OFsDYZwCVsRtM0EpITSLqQRDvHdvy2/2/xcJK6wkJYmwyLEDazPW07SReScLBz4IHwByTpC9FEJPELm9iXuY/vMr7DTmfHjNAZ+La7/sK/EMI6pMSSaHLHLhxjS+oWACb3mUxvr953tkGTSdUTuHz56lJWVrtPnL391cIxrq7g6anqEUuVMdEGybdeNKnkvGQSkhMAGBc47vbmzC0oULOOVJeYzMlRyf9W6XSq9rC3txoa6ecHBoPUIhatnlzcFU0mOS+ZNSfWYNJMDO0+lNEBDey2aTKpRJ+SAr/8ohL/tTp0UGfx1TUEXF1rj3uvqlI1A6rrBhQW1l2PWKdTRekDAiAwUE1fJb8KRCsjiV80iVN5p/j8xOeYNBODDYMZ22ts/d02NU2d0R87pmYZLy29+pqLi0rMfn7qTL1Ll9s7SzcaIT8fLl6ErCzIzITz52uPx3dwULOZ9+mjFlfXW/8cIZoZSfzC6n6++DOfn/wcs2bmLr+7uCfgnhsn/YoKlex/+EEl5Gre3tC3rzoL79bNenV6q6rUQSA1Vf26uHDh6mt2duogEBoqBwHRokniF1b1Y86PfHXqK8yamSF+QxgTMKbupF9QAAcOQFLS1Rq+7u4QHg4REeqs3hYDuy5fhtOn4eRJSEu7+mvAzk796oiIgJAQ9ctAiBZCEr+wCk3T2J+139J75+7udzO65+jrk35BAXz7LRw/fjWp9ugBAweqs+rmVBSltBROnYITJ2ofBJyc1K+RyEgVu4w8Fs2cJH7R6Myamc2/bOZAtirDHdcrjsF+g2uvlJ9/NeFrmjqDjoiAQYOsU46wsZWWqgPAsWOqaaha+/bqABAZqSaoFaIZksQvGlWVqYqvTn3FyYsnsdfZc3/I/YT5hF1doa6EHxWlZpvo0MF2gd+JvDx1ADh+XE1DVa17d3UACA1VcwgK0UxI4heNprCskNUnVnPhygWc7J2YFTaLnh16qhdbY8K/lqZBerq6TvHzz2pGEFDdQfv0UX9vQIBMICtsThK/aBS/FPzCFye/oMxYRkeXjswKm4WPm0/bSPh1qaxUF4SPHVPXA6q1a6eatKKiVE8lIWxAEr+4I2bNzJ6ze9iVvgsNjSCvIKaETMG5uKRtJvy6FBWp/ZCUVHvwma+v2idhYdI1VDQpSfzithWWFfLVqa/IKM5Ah44R/iMY5tYX3Z498OOPkvCvpWlqkFj1oLTqbqv29hAUpPZTYGDz6skkWiVJ/OKWaZrGkfNH2Jy6mUpTJe0c2zHJ524Cj2fXTvjR0Srht29v65Cbn6oqSE5WvwJSU6+WjnB1VV1D+/YFf3+5HiCsQhK/uCUFZQVsSNnALwW/ABCq78qETBdck8+oFSTh37rLl682BdUcrezmpgaHhYaq8QFyEBCNRBK/aJAqUxV7MvawN2MvJrMR5yvlTMj1JDyzRs+V6GgYMkQS/u3SNMjNVeMDTpxQF8arubmpXwEhIeogIM1B4g5I4hf10jSNkxdPsvXMVopKC+DiRaIu2jMmzwN3HMHREWJi1MCrdu1sHW7roWmq3HT1QaDmRWEnJ3UtICgIeveWC8PilkniF3XSNI3k/GR2pu0kpzATzp2jy4XL3Fvene54qgqZsbFqcXGxdbitm6apYnEnT6rrArm5V1/T6dRAsV691NK1qzQJiZuSxC9qMWtmkvOS2XP2W85ln4Jz5/C4WMxwc3ei6YqdT2dVRyciQp3ti6ZXWKgKxyUnw9mztSehcXFRFUQDAtSBoH17qR0kriOJXwBQYazg6IWjHPhlF4Vnk+HCBdzLTAylO/11vuiD+6qze39/SSTNSXm5mqQmNVXdFhbWft3DQ/0i6N5dXRvw8ZF/PyGJvy3TNI2M4gyOpe3nxM/fUpGTDZeK6aA5E4uB/u5BOET1gwED5IJtS1FYePUgkJam5h6uydlZTS/ZrZsaQObrK9dm2iBJ/G2MpmnkXsnh55R9HDu1m8KcdEthMX/aM8i+B0F9hmAX3U/qyrR0mqYKyJ09CxkZ6rZmEblq7dpdPQj4+kLnzuo5+WXQaknibwOMZiNZWSc5nbyPUxlHKMg9axk16oETEfa+RAbchXfEINVTRCYbb72KiyE7Wy3nzqmpJsvLr1/P2VnVEvLxqX3r7i4HhFZAEn8rZDRWkpN5ivTUw6RlHONsbgpV5SWW191wINjZQGhALD1D78YusLck+7ZK01RX0XPnri65udc3EVVzdlbzDHTsqEpw1LyVXwkthiT+Fs5cXkZe1mnOZZwg+3wy2Xlp5Fw6h8lkrLWej4MngT4h9AmIwdAnFjvfbvKfVNRN06CkRB0ALl6sfVvXr4Nqer26FuThoRZPz6v3qxdnZ/neNQOS+Js7TUMrKaGsMIeivGzyctPJy8/gYmE2eZdyKKgoxETtf0Id4O3SCUOnAHr2iKRn8CDcu/WU9npxZ6oPCAUF6iLytbclJTffhl6vRiFfu7i7X73v6qq6pTo7q1+icqBodFZN/Js2beKpp57CZDKxYMECFi1aVOv1iooK5syZw+HDh/Hy8mL16tX4+/tbK5zmQ9OgogKtpISKSwWUXM6n5FI+JSWFlFwppKS0mEtX8ii6kk9xWSHF5lKqMNe9LTs72rt2xLdjD7p16U03v1C6+ofj1E564YgmVlGhriFcunT1tnqpflw9OU1D6XQq+VcfCGou1c85Ol5dHBxufF+vl4PIr6yW+E0mE0FBQWzduhWDwUBMTAwrV66kb9++lnXee+89jh8/zvvvv8+qVav46quvWL16tTXCqZumqcVsrn8xmTBXVWKsLMdYVaGWX+9XVVVgNFZYnq+qqqCysoyKihIqKkqpqCylorKMiqoyKqrKqTCqpUyropSq687W66TX4+zsjqdrB7w8uuDt5Ucnb386dQ3Ay8cfRweZ1k+0EJWV6pfBzZbycrVUl65uDDrd1YNA9YFAr1d1j271vr29+gVdc7nd53S662/req4R6Rt1azUcPHiQwMBAAgICAJg1axYJCQm1En9CQgIvvvgiANOmTeOJJ55A0zR0jfBHlhXn8/UHv8dsNmPWTJg189Wl5nNoN11Mv942Knt7cHDCydEVNyd33Jw9cHPxwM3VE1dXTzw8ffDs6Iunly+eHj446yW5i1ag+gy8oXMzmM0q+ZeVXT0YVC/Vz1VWqqWqqu771Y+NRrWtxjyYNIXZs1VNpkZktcSfnZ2Nn5+f5bHBYODAgQM3XEev1+Pp6Ul+fj6dOnW648938fTigWf/c8fbEULYkJ2datKRelCNympX++pqQbr2TL4h6wghhGhcVkv8BoOBzMxMy+OsrCx8fX1vuI7RaKS4uJiOHTtaKyQhhBBYMfHHxMSQkpJCWloalZWVrFq1ivj4+FrrxMfH8/HHHwOwdu1aRo0aJWf8QghhZVZL/Hq9nsWLFxMXF0dISAgzZswgNDSU559/nsTERAAefvhh8vPzCQwM5J///CevvfZao33+pk2bCA4OJjAwsFG3e6syMzMZOXIkISEhhIaG8q9//QuAgoIC7rnnHnr37s0999xD4bVVFZuIyWQiOjqaiRMnApCWlkZsbCy9e/dm5syZVN5q97tGUFRUxLRp0+jTpw8hISHs37+/Weyv//3f/yU0NJSwsDAeeOABysvLbba/5s+fj4+PD2FhYZbnbrSPNE3jySefJDAwkIiICI4cOdKkcT3zzDP06dOHiIgI7r//foqKiiyvvfrqqwQGBhIcHMzmzZubNK5qb775Jjqdjry8PMD2+wvgnXfeITg4mNDQUJ599lnL8422v7RWyGg0agEBAVpqaqpWUVGhRUREaCdOnLBJLOfOndMOHz6saZqmXbp0Sevdu7d24sQJ7ZlnntFeffVVTdM07dVXX9WeffZZm8T3j3/8Q3vggQe0CRMmaJqmadOnT9dWrlypaZqmLVy4UHvvvfeaPKY5c+ZoH374oaZpmlZRUaEVFhbafH9lZWVp/v7+WmlpqaZpaj8tW7bMZvtr9+7d2uHDh7XQ0FDLczfaR+vXr9fGjRunmc1mbf/+/drAgQObNK7NmzdrVVVVmqZp2rPPPmuJ68SJE1pERIRWXl6unTlzRgsICNCMRmOTxaVpmpaRkaGNHTtW6969u3bx4kVN02y/v3bs2KGNHj1aKy8v1zRN03JycjRNa9z91SoT/759+7SxY8daHv/tb3/T/va3v9kwoqvi4+O1LVu2aEFBQdq5c+c0TVMHh6CgoCaPJTMzUxs1apS2fft2bcKECZrZbNa8vLws/0mv3Y9Nobi4WPP399fMZnOt5229v7KysjSDwaDl5+drVVVV2oQJE7RNmzbZdH+lpaXVShg32kePPPKI9tlnn9W5XlPEVdOXX36p/eY3v9E07fr/l2PHjtX27dvXpHFNnTpVS0pK0nr06GFJ/LbeX9OnT9e2bt163XqNub9aTXq7mwAAB6dJREFU5Rj+urqSZmdn2zAiJT09naNHjxIbG0tOTg5du3YFoGvXruTWnE6vifz+97/n73//O3a/lnLIz8+nffv26PWql68t9tuZM2fw9vZm3rx5REdHs2DBAkpKSmy+v7p168Yf//hHunfvTteuXfH09KR///4231813WgfNaf/D0uXLmX8+PHNIq7ExES6detGZGRkredtHdfp06fZs2cPsbGxDB8+nB9++KHR42qViV9rht1Er1y5wtSpU3nrrbfw8PCwaSwA69atw8fHh/79+1ueaw77zWg0cuTIER577DGOHj2Km5ubTa/RVCssLCQhIYG0tDTOnTtHSUkJGzduvG49W3/P6tIc/l0BXnnlFfR6PbNnzwZsG1dpaSmvvPIKL7300nWv2Xp/GY1GCgsL+f7773njjTeYMWMGmmqdabS4WmXib0hX0qZUVVXF1KlTmT17NlOmTAGgc+fOnD9/HoDz58/j4+PTpDHt3buXxMRE/P39mTVrFjt27OD3v/89RUVFGI2qsqct9pvBYMBgMBAbGwuoEd1Hjhyx+f7atm0bPXv2xNvbGwcHB6ZMmcK+fftsvr9qutE+ag7/Hz7++GPWrVvHihUrLMnKlnGlpqaSlpZGZGQk/v7+ZGVl0a9fPy5cuGDz/WUwGJgyZQo6nY6BAwdiZ2dHXl5eo8bVKhN/Q7qSNhVN03j44YcJCQnh6aeftjxfsyvrxx9/zKRJk5o0rldffZWsrCzS09NZtWoVo0aNYsWKFYwcOZK1a9faLK4uXbrg5+dHcnIyANu3b6dv374231/du3fn+++/p7S0FE3TLHHZen/VdKN9FB8fzyeffIKmaXz//fd4enpamoSawqZNm3j99ddJTEzE1dW1VryrVq2ioqKCtLQ0UlJSGDhwYJPEFB4eTm5uLunp6aSnp2MwGDhy5AhdunSx+f6aPHkyO3bsAFSzT2VlJZ06dWrc/XVbVwZagPXr12u9e/fWAgICtJdfftlmcezZs0cDtPDwcC0yMlKLjIzU1q9fr+Xl5WmjRo3SAgMDtVGjRmn5+fk2i3Hnzp2WXj2pqalaTEyM1qtXL23atGmWngVN6ejRo1r//v218PBwbdKkSVpBQUGz2F/PP/+8FhwcrIWGhmr/9V//pZWXl9tsf82aNUvr0qWLptfrtW7dumkfffTRDfeR2WzWHn/8cS0gIEALCwvTfvjhhyaNq1evXprBYLB8/xcuXGhZ/+WXX9YCAgK0oKAgbcOGDU0aV001L+7aen9VVFRos2fP1kJDQ7Xo6Ght+/btlvUba39JPX4hhGhjWmVTjxBCiBuTxC+EEG2MJH4hhGhjJPELIUQbI4lfCCHaGKvNwCWEreTn5zN69GgALly4gL29Pd7e3gC4urqyb98+q3xu//792b9/P46Ojpbnli9fzqFDh1i8eLFVPlOI2yGJX7Q6Xl5eJCUlAfDiiy/i7u7OH//4R6t+Znp6Ot26dauV9IVorqSpR7Qp7u7uAOzatYvhw4czY8YMgoKCWLRoEStWrGDgwIGEh4eTmpoKwMWLF5k6dSoxMTHExMSwd+/eOre7ceNGxo0bB8CyZcsICgpi+PDhtdb/5ptviI2NJTo6mjFjxpCTk4PZbKZ3795cvHgRALPZTGBgIHl5eXz++eeEhYURGRnJsGHDrLlbRFtz52PPhGi+XnjhBe2NN96wPHZzc9M0TY1W9vT01M6dO6eVl5drvr6+2v/f3t27pBbHYQB/sMHJTWgIh2yJ0IggxN5oCEeFICVoKLCppJc/IA7R0Fa0ZARtTS2NThmcIQIHIWgIGsKTQWDRy1AGPXc49ONKL3BvXu5wns/kkZ8edXiG30+e7/LyMklyY2OD8/PzJMmJiQnatk2SvLy8ZGdn56f3SSaTvLi4YLVaZSgU4s3NDV9eXtjf38/Z2VmS5O3tramb3tnZ4dLSEknSsiyur6+TdLvrx8bGSJKRSISO45Ak7+7umvejiOdpq0c8q6+vz3SwdHR0IJFIAHB7XIrFIgC3nO3s7My85uHhAY+PjwgEAua5er0Ox3EQDodxcHCAkZERc6aQyWRwfn4OwC3VymQyuL6+Rr1eR3t7OwB3ClMqlcLCwgJ2d3cxPT0NABgYGMDU1BTS6bQp9xNpBm31iGf5/X7z2OfzmWufz2caN9/e3nB8fIxyuYxyuYyrq6uG0AcA27YxODhorr+qys3lcpibm8Pp6Sm2t7fx/PwMAAiFQmhtbcXh4SFOTk5MX30+n8fq6ioqlQp6enpQq9Wa9+XF0xT8It9IJBIN/8h5PzT+XaFQMGEdi8VwdHSEWq2G19dX7O/vm3X39/doa2sDANOi+S6bzWJychLpdBotLS0A3OrgWCyGlZUVBIPBhkpekZ9Q8It8Y3NzE6VSCd3d3ejq6kI+n/+w5v2gGHAnX1mWhXg8jtHRUfT29pp1lmVhfHwcQ0NDCAaDDe+RTCbx9PRktnkAd0h5NBpFJBLB8PDwh0lRIn9L7ZwiP+A4DmZmZj6dxvUnSqUSFhcXYdt2kz6ZyNcU/CL/2draGra2trC3t9dwViDyryj4RUQ8Rnv8IiIeo+AXEfEYBb+IiMco+EVEPEbBLyLiMQp+ERGP+QUTbS7Azc27TAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the data on three separate curves for S(t), I(t) and R(t)\n",
    "fig = plt.figure(facecolor='w')\n",
    "ax = fig.add_subplot(111, axisbelow=True)\n",
    "ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')\n",
    "ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')\n",
    "ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')\n",
    "ax.set_xlabel('Time /days')\n",
    "ax.set_ylabel('Number (1000s)')\n",
    "ax.set_ylim(0,1.2)\n",
    "ax.yaxis.set_tick_params(length=0)\n",
    "ax.xaxis.set_tick_params(length=0)\n",
    "ax.grid(b=True, which='major', c='w', lw=2, ls='-')\n",
    "legend = ax.legend()\n",
    "legend.get_frame().set_alpha(0.5)\n",
    "for spine in ('top', 'right', 'bottom', 'left'):\n",
    "    ax.spines[spine].set_visible(False)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
